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ABSTRACT 

In a disc galaxy the distribution of azimuthal components of velocity is very skew. In the past 
this skewness has been modelled by superposed Gaussians. We use dynamical arguments to 
derive an analytic formula that can be fitted to observed velocity distributions, and validate it 
by fits to the velocities derived from a dynamically rigorous model, and to a sample of local 
stars with accurate space velocities. Our formula is much easier to use than a full distribu- 
tion function. It has fewer parameters than a multi-Gaussian fit, and the best-fitting model 
parameters give insight into the underlying disc dynamics. In particular, once the azimuthal 
velocities of a sample have been successfully fitted, the apparatus provides a prediction for 
the corresponding distribution of radial velocities Vr, An effective formula like ours is invalu- 
able when fitting to data for stars at some distance from the Sun because it enables one to 
make proper allowance for the errors in distance and proper motion when determining the 
underlying disc kinematics. The derivation of our formula elucidates the way the horizontal 
and vertical motions are closely intertwined, and makes it evident that no stellar population 
can have a scale height and vertical velocity dispersions that are simultaneously independent 
of radius. We show that the oscillation of a star perpendicular to the Galactic plane modifies 
the effective potential in which the star moves radially in such a way that the more vertical 
energy a star has, the larger is the mean radius of its orbit. 

Key words: galaxies: structure - kinematics and dynamics - Galaxy: disc - solar neighbour- 
hood - kinematics and dynamics - stars: kinematics 



1 INTRODUCTION 

Currently considerable effort is being invested in surveys of the 
solar neighbourhood. Fifteen years ago the study of nearby stars 
was revived by the Hipparcos mission, which pioneered space as- 
trometry. Hipparcos put ground-based astrometry onto a more se- 
cure foundation, so now useful proper motions are available for 
tens of millions of stars. In the last decade the proper motions have 
been complemented by photometric surveys, both in the infrared 
and to fainter magnitudes at optical wavelengths. Finally, the RAVE 
(Steinmetz et al. 2005) and SEGUE (see York et al. 2000; Yanny 
et al. 2009) surveys have measured nearly a million line-of-sight 
velocities. 

As a result of these major observational programmes, it is be- 
coming possible to determine the velocity distribution within the 
disc of our Galaxy, not only at the location of the Sun, but also at 
significant distances, especially at higher Galactic latitudes. Nat- 
urally one wants to quantify the velocity distribution observed at 
some location x in the Galaxy in an efficient way. Conventionally 
one does this by imagining that the density of stars in velocity space 
forms a "velocity ellipsoid" - a triaxial ellipsoidal region of over- 



E-mail:rasch@mpa-garching.mpg.del 



density in velocity space. If the Galaxy were axisymmetric (which 
is a reasonable first approximation), we would expect that in the 
Galactic plane the principal axes of the the velocity ellipsoid would 
be aligned with the coordinate directions of cylindrical polar coor- 
dinates, (R, z, </>)• As one moves above the plane, two of the princi- 
pal axes of the velocity ellipsoid are expected to tip slightly with 
respect to the R and z directions. Let the components of velocity 
parallel to these principal axes be denoted vi and v 2 , where V\ — > v R 
and \>2 — > v ? as z — > 0. The third axis is expected to remain aligned 
with the (j) direction. 

The distributions of the Vi and v 2 components of velocity are 
expected to be roughly Gaussian with vanishing means and to be 
to good approximation aligned with the Galactic polar coordinates 
(minor vertex deviations as found in the solar neighbourhood by 
Dehnen 1998, will not be discussed here). Consequently, they can 
be characterised by their standard deviations a\ and <x 2 . The distri- 
bution of the V0 components peaks at a value of that is slightly 
smaller than the circular speed v c . However, it is not at all well 
modelled by a Gaussian, because it is very skew, with many more 
stars at v$ = v c — v than at v c + v, causing the population to have 
non-zero asymmetric drift. Notwithstanding this skewness that was 
already known to Gustav Stromberg (Stromberg 1927), distribu- 
tions have traditionally been characterised by a mean and a stan- 
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dard deviation. Since a single Gaussian fits the data very poorly, 
the observed distribution is frequently modelled by a superposition 
of two Gaussians: then the overall distribution is characterised by 
two means, two dispersions and the ratio of the numbers of stars 
accounted for by each Gaussian, a total of five shape parameters. 

The purpose of this note is to introduce a new representation of 
distributions that is more effective in the sense that it fits typical 
data more accurately with fewer and physically more meaningful 
parameters. Moreover, the new representation has a dynamical ba- 
sis, so it is able to connect the skewness of the distributions to 
the standard deviations in v\ and v 2 . With the new representation, 
a single free shape parameter suffices to describe the distribution 
in for the whole population of stars in the solar cylinder, and 
two parameters are sufficient to fit the distribution in of stars that 
have a given distance from the plane. The formula predicts in a nat- 
ural way both the magnitude of the asymmetric drift and the offset 
of the modal azimuthal velocity from the circular velocity, neither 
of which is achieved by multiple Gaussians. 

The paper is organised as follows. In Section 2 we derive an 
approximation to the distribution in an annulus in the Galactic 
disc that can be used as standard for extragalactic measurements. 
As most samples of the Galaxy are centred on certain Galactic alti- 
tudes |z|, in Section 3 we use the adiabatic approximation (Binney 
2010) to take into account the vertical motions of stars. In 3.1 we 
derive a formula that accounts for the variation in the v<» distribu- 
tion with \z\ and test it against the velocity distributions of more 
elaborate models. In Section 3.2 we derive a formula for the way 
in which the in-plane motion of a star depends on the extent of its 
excursions perpendicular to the plane. The outcome is a small cor- 
rection to the distribution derived in Section 3.1. In Section 3.3 
we give formulae from which the distributions of v R and v z follow 
once the distribution in has been fitted. In Section 4.1 we show 
that our formulae provide good fits to the disc model of Binney & 
McMillan (2011; hereafter BM11), which has a rigorous dynam- 
ical basis. In Section 4.2 we demonstrate the practical application 
of the formula by fitting data from the Geneva-Copenhagen Survey. 
Section 5 sums up and looks to the future. 



2 VELOCITY DISTRIBUTION AS A 2D PROBLEM 

There are three reasons for the asymmetry of the distribution of 
components of nearby stars: stars at low v$ are approaching apoc- 
entre, so they have guiding-centre radii R g smaller than the solar 
radius, R . As one moves inwards through the disc, not only does 
the density of stars increase rapidly on account of the exponential 
increase in the surface density l.(R) oc expi-R/Rj), but the random 
velocities of stars also increase, so a greater fraction of all stars 
are on eccentric orbits that carry them far from their guiding-centre 
radius R g . Moreover, the effective potential in which a star oscil- 
lates around R g rises much more steeply at R < R g than it does at 
R > R g , so stars spend more time beyond R g than they do interior 
to it. In fact, as a population of stars heats up over its lifetime, the 
asymmetry of the effective potential causes the population to ex- 
pand spatially, and by conservation of angular momentum its mean 
rotation rate diminishes. For all these reasons, there are many more 
visitors reaching Rq with guiding centres at Rq - A than at Rq + A. 
A functional form for n(v^) that is successful in fitting observed 
distributions will reflect these facts. 

Following Shu (1969) we decompose the energy of a disc star 



into three parts. If the star were on a circular orbit with angular 
momentum L z , it would have energy 

E c (L z ) = <b eB {R e ,L z ), (1) 
where with the Galactic potential in the plane 0(R) 

D eff (*,4)= (2) 

and the guiding-centre radius R g (L z ) solves the equation R g v c (R g ) = 
L z , with v c (R) the circular speed. In addition to this energy, the star 
has two smaller energies, namely the energy E z of vertical motion 
and the energy E R of random motion within the plane. We postpone 
discussion of E z and focus for now on stars with E z = 0, which 
move in the plane. We have 



E R = k V R + ®eff (R, L z ) - <t> eS (R g ,L z ) 
= ^ 2 R + A^(R,L Z ) 

where 

A(D cff (fl,L z ) = ® ee (R,L z ) - *rff(fi E ,L z ). 



(3) 



(4) 



Suppose that the disc's distribution function (df) is (Shu 1969) 



f(E R ,L z ). 



_ p -Er/<t 2 



(5) 



where F(R g ) is a function that determines the surface density of 
the young disc and cr(R g ) is a function that determines how the ra- 
dial velocity dispersion varies with R. On account of the tendency 
noted above for a population to expand radially as it heats up, if 
F(R g ) is the same for both cool and hot populations, the hotter 
populations will have slightly larger radial scale-lengths than the 
cool ones. Note that cr gives the intrinsic dispersion of stars at their 
guiding centre radius R g . The dispersion <v|> 1/2 actually measured 
at some radius R will have contributions from all populations that 
reach this radius and turns out to be -10 per cent higher than the 
value of cr for the stars that have R g = R on account of the presence 
of stars that have guiding centres at R g < R. 

Schonrich & Binney (2009) show that the probability per unit 
area that a star with angular momentum L z will be found at R is 



P(R\L Z ) = JL exp 
<tk 



A<t> eff (R,L z ) 



(6) 



where K(R g ) is chosen such that 1 = In fdRRP. 

Let n(v$, R) dv^ be the number per unit area of stars at R with 
in (v^,v^ + dv ). Then 

n(v„ R) dv, = N(L Z ) dE z P{R\L Z ) = N(L z )P(R\L z )Rdv^ (7) 

where N(L,)dE z is the number of stars in the disc with L z in 
(L z , L z + dL z ). In a cold disc, the number of stars with angular mo- 
menta in (L Z ,L Z + dL z ), is simply the mass in the corresponding 
annulus, 2nl,(R g )R g dR g . Hence under the neglect of the mild radial 
expansion noted above of a population as its dispersion increases, 
an exponential disc with surface density 

Z(R) = T <) e <R - R o ySi , (8) 

where S is the local surface density and R d is the radial scale- 
length of the disc, has 



N(L Z ) : 



2n'S. Q R g 



-(Rg-R )/Kd 



where we have assumed a constant circular speed, so dL z 
Combining equations (6), (7) and (9), we have 



n(\^,R) 



v c cr 



A<t> eff (R,L z ) 



(9) 
v Q dR g . 

(10) 
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Our assumption of constant v c allows us to evaluate K(R„), 
because then 



-A<D cff (/?,L z ) = \_L 2 z (Rf 2 - R- 2 ) + v 2 c ln(R g /R), 

so the normalisation condition 1 = 2n J dR RP reads 



(11) 



— = f 



dR exp 



AO cff (tf,L z ) 



dR sxp[c{2\nR g /R+l-R g 2 /R 2 )] 
= g(c)R g , 



where 



c(R g ) = 



20-HK) 



and 



8(c) ^ 



So 



e c (c 



1)1 



2<*-i/2> 



(12) 



(13) 



(14) 



Vcg(c) 
So 

Vcg(c) 



exp 



exp 



R e -R t 



c\2\n— + 1 
R 



R 2 



Rn — R(, 



(15) 



As here we are aiming at velocity distributions and not stellar den- 
sities at a certain position, we will henceforth use the normalised 
velocity distribution at a fixed radius R 



N 

n(v,p\R) = ——exp 
g(c) 



H21n-^ + 1 



Rf 
~R 2 ~ 



R B -R 



Rd 



(16) 



where N normalises the integral of n in v<j to unity. 

Note that on the right side of eq. (16) the dependence on is 
carried by the instances of R g = Rv,p/v c and by c(R g ). 

From equation (6) we expect K (which has the units of a fre- 
quency) to be constant when <r « v„ and indeed from an asymp- 
totic expansion of equation (12) for large c we obtain g(c) oc <r'' 2 oc 
cr and K = v c /(2n 3/2 R g ). Fig. 1 shows that g(c) satisfies this ex- 
pectation throughout the entire parameter range of interest. For 
(c > 25), the direct computation of g(c) becomes impractical so 
apart from the advantage of a numerically less costly formula a 
reasonable approximation must be found. The dashed green line in 
Fig. 1 demonstrates that for c > 2 this is achieved to high precision 
by 1 



g(c) ■■ . x , 

A/ 2(c- 0.913) 

For cr we adopt the radial dependence 

<r(R g ) = o- Qe -^- R °^ . 



(17) 



(18) 



Above we have restricted ourselves to stars with E- = 0. How- 
ever, to the extent that the motion in R of a star is unaffected by its 
motion perpendicular to the plane, the distribution we have derived 
will apply to the population formed by all stars that now lie in the 



1 Alternatively g(c) can be stringently approximated using eq. (8.327) of 
Gradshteyn & Ryzhik (1980), but our term appears to be the best com- 
promise of simplicity and accuracy throughout the interesting part of the 
parameter range. 




0.00001 



Figure 1. Upper panel: The function g(c) defined by eq. (14) with solid red 
line and the approximation of eq. (17) shown by a dashed green line. Lower 
panel: the relative difference between the approximation and the underlying 
formula. 



solar cylinder (the region restricted in radius to R = Rq but unre- 
stricted in z). From our formulae we have that the shape of this ve- 
locity distribution is controlled by four parameters: the galactocen- 
tric radius of the measurement Rq, the scale-length of the young 
disc, the local velocity dispersion cr , and the scale-length R^ on 
which the velocity dispersion varies. The first two parameters are 
generally well-known and for fits of the solar neighbourhood can 
be set to R = 8 kpc and R A = 2.5 kpc. The value of R a is less clear 
and will be discussed below. 

The dependence of n(v,p\R ) on <x and R a is shown in Fig. 2. 
In the upper panel we hold the dispersion scale-length constant at 
7.5 kpc and show the velocity distributions for local dispersion val- 
ues of <r = 20, 25, 30, 40, 50km s~' . As cr increases, the distribu- 
tion becomes wider and the low-velocity tail rises much faster than 
does the high-velocity tail, corresponding to an increasing asym- 
metric drift. Simultaneously, the peak slowly shifts to lower ve- 
locities. The lower panel shows the velocity distributions for fixed 
<r = 30kms _1 , but different scale-lengths of the velocity disper- 
sion, R^ = 10, 7.5, 5, 4 kpc. Smaller values of R a imply higher dis- 
persion in the inner regions of the Galaxy and lower dispersion in 
the outskirts. Thus for small R a stars from the inner disc can more 
easily reach the R compared to their counterparts with larger R g 
so the azimuthal velocity distribution is more skewed and develops 
a strong low-velocity tail. In the most extreme case, R^ = 4 kpc, 
the model breaks down, as (v^) 1/2 approaches v c in the inner re- 
gions. Notice that as R a falls, the mode of the velocity distribu- 
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Figure 2. The behaviour of eq. (16) for different local dispersions (upper 
panel) and different scale-lengths on which the dispersion varies (lower 
panel). The radius of observation is R = 8 kpc at a young-disc scale-length 
Ri = 2.5 kpc. In the upper panel we set the scale-length of the velocity 
dispersion to R^ = 3R<i = 7.5 kpc and show results for local dispersions 
o"o = 20, 25, 30, 40, 50km s _1 , while we use ctq = 30kms~' in the lower 
panel. 



tion shifts to higher even as the tail at high becomes weaker. 
This effect reflects the fact that stars with large random velocities 
spread their contributions to the velocity distribution over many 
radii, while stars with small random velocities localise their contri- 
butions. Consequently, when R a is small, the width in R over which 
stars contribute to narrows strongly as R g increases, which puts 
stars with R g = Rq strongly in control of the mode of the local 
distribution. 



3 THE VELOCITY DISTRIBUTION AS A FUNCTION OF 
DISTANCE FROM THE PLANE 

In the previous section we set E- = to obtain results that are ap- 
proximately valid for the distribution of stars in regardless of 
their distance from the plane. In practice we can determine the ve- 
locity distributions of the stars that lie in various more-or-less nar- 
row ranges in z. We must now consider how these distributions will 
vary with z and differ from the aggregate distribution determined 
above. In the following we do so using a number of quite crude 
physical approximations. These approximations give useful physi- 
cal insight into why the distribution in varies with \z\ as it does, 



but ultimately the value of our final fitting formula does not depend 
on the correctness of the arguments used to motivate it. 



3.1 Weights of different populations as functions of z 

The key idea is that within the solar cylinder there coexist many 
populations, one for each value of R g . Indeed the chemical compo- 
sitions and ages of stars vary systematically with R g (see e.g. Luck 
& Lambert 201 1; Bensby et al. 201 1; Schonrich, Binney & Dehnen 
2010, for observations and a short discussion of consequences for 
studies of kinematics). Moreover, the smaller a population's value 
of R g , the larger will be its mean value of E z and therefore the larger 
will be its vertical scale-height h at radius R; here h(R g ,R) is the 
distance that at radius R provides the best fit to the vertical density 
profile of the population through 

n R AR,z)oc ye- ulh . (19) 
s h 

Note that with this formula we are not asserting that the real vertical 
density profile is exponential, but simply identifying the character- 
istic vertical extent of the population. We now investigate how the 
vertical extent of the population increases with R because the verti- 
cal restoring force, which scales like l.(R), decreases outwards. 

BM11 show that a very good approximation to the vertical 
dynamics of a population of stars can be obtained by assuming that 
the vertical action 



J- 



dzv z 



(20) 



of the population's stars is adiabatically invariant as the stars oscil- 
late in radius. We use this adiabatic approximation (aa) to estimate 
the ratio h(R g ,R)/h(R g ,R g ). 

J z can be evaluated analytically only for a vertical force K z that 
is proportional to z" and to the local surface density of the disc E, 
so we assume these dependencies in order to gain an analytic model 
- in reality K z has a much more complex dependence on z, which 
does not yield an analytic expression for E Z (J Z ). Then the vertical 
action is 



J, = 



2 3/2 fJm 

dz yjE- - KLz\ 

" Jo 



where A; is a constant and 



(EJkZ) 



l la 



(21) 



(22) 



is the height at which the radical vanishes. In terms of the variable 
9 = z/z m we have 



J- 



2 3/2 

n 



f AO Vl -0". 
Jo 



(23) 



With the aa, it now follows that 

z m oc 5T'/< 2+ ">. (24) 
The scale-height h(R g ,R) of the population will scale like z m , so 
h(R g ,R) I Z(R) \- [/>2+a) I R-R„ 



h(R g ,R g ) 



\m g )l 



■ exp 



(2 + a)R d 



(25) 



Photometry of edge-on spiral galaxies shows that the overall scale- 
height varies very little with radius (van der Kruit & Searle 1982). 
Since the velocity dispersion at radius R will be dominated by the 
population that has R g = R, we assume that 



h = h(R,R) 



(26) 
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Figure 3. The relative contributions given by eq. (30) of the populations 
with guiding-centre radius R g measured at galactocentric radius R to the 
velocity distribution at various values of |z|. We use a = 1.5, /io = 300pc 
and Rj = 2.5 kpc. 



is independent of R. Hence 



h(R g ,R) 
h(R, R) 



exp 



(2 + a)R A 



(27) 



For stars that make only small-amplitude vertical oscillations, 
K, oc z so a = 2. If the amplitude of a star's oscillations signifi- 
cantly exceeds the local disc scale-height (but is none the less small 
compared to R), a better approximation is that the disc is razor thin, 
so K z ^ constant and a = 1 . For amplitudes not small compared to 
the disc's scale-length, a yet smaller value of a is appropriate. In the 
fits described below we assumed that a decreases with z according 
to 

'2- 1.5z/1.5kpc forz<1.5kpc 
, 0.5 otherwise. 



a{z) - 



(28) 



These choices are educated guesses that produce useful results. 

Let the ratio of the contributions to the distribution in at 
height z of the population with guiding-centre radius R g and the 
local population be given by the factor f(z,R g - R)- Then from 
equation (15) the distribution of at altitude z above the plane is 



N 

n(y^\R,z)= —-exp 
8(c) 



-(R,-R y 




AO cff 




exp 









f(z,R g -R), (29) 



where N is a new normalisation constant. 

By equation (19), the factor / will scale as 



f(z, R g -R)- 



n R (z) 
h(R, R) 



■ exp 



h(R g ,R) 

R. - A' 

= exp 



Mil- 

«0 



h(R,R) 
h(R g ,R) 



(30) 



, exp < — 1 - exp 
(2 + a)R d ) V \ho X ' 



R g -R 
(2 + a)R A 



where in the second step we have used the definition (26). Fig. 3 
shows, for five values of z, how / varies with R g at R = R . At z = 0, 
/ is an increasing function of R g because at large R g stars typically 
have small E- and therefore are more numerous in the plane than in 
the solar cylinder as a whole. At z = 1 kpc, by contrast, / falls with 
increasing R g because high above the plane a sample is richer in 
stars with small R E , and therefore large E z , than is the solar cylin- 
der as a whole. The formula attains a maximum, where the local 
scale-height of a population equals the altitude. Since by definition 



f(z, R g - Rq) gives the enhancement of stars with guiding-centre ra- 
dius R E relative to local stars, all lines intersect at (/ = 1, R g = i?o) 
irrespective of the chosen altitude. 



3.2 Impact of the AA on the radial motion 

The aa predicts that any star's value of E : varies as it moves radi- 
ally. Since the star is moving in a time-independent potential, its 
total energy is constant. It follows that changes in E- must be com- 
pensated by changes in its energy of motion parallel to the plane. 
Since L. is a true invariant, the energy O e fj(,R g , L ; ) associated with 
the underlying circular orbit is unchanged, so changes in E, must 
be compensated by changes in the radial energy E R . Hence (3) be- 
comes 



\ v\ + A<D e ff(^, L z ) + AE Z = constant, 
where 

AE z (J z ,R,R g ) = E Z (J Z ,R) - E z (J z ,R g ) 



(31) 



(32) 



takes into account the decrease in a star's vertical energy as it moves 
outwards: by conservation of energy, this energy must be trans- 
ferred to the radial motion. One way of thinking about this energy 
transfer is to imagine that the star moves at constant energy in an 
effective potential 



AO ad (tf,L,) = AO cff (,R, L z ) + AE z (J z ,R,R g ), 



(33) 



which might be called the "adiabatic potential" since it is the ef- 
fective potential for radial motion that follows from adiabatic in- 
variance of the vertical motion. At R > R, AO a d increases with 
R less rapidly than the standard effective potential AO c ff, so stars 
can reach larger radii than they could if E z were constant. BM11 
simulated this effect by simply increasing L z to L z + J z . 

With the power-law forms of the vertical potential that we in- 
troduced above, we can obtain the dependence of E z on X. From 
equations (22) and (23) we find 



2/2+ff 



AE z (J.,R,R g ) = E Z (J Z ,R) 



1 - exp 



2(R - R g ) 
(2 + a)R A 



(34) 



(35) 



Our plan is to use AE Z to modify our expression (10) for ^(v^l^, z), 
which has no dependence on v ; and therefore does not specify 
a value of J- or E z . Therefore in equation (35) we now replace 
E Z (J Z ,R) by an estimate E z of the typical vertical energy of the stars 
that are encountered at (R,z) with angular momentum L, = Rv$. 
Equation (47) in the Appendix is an expression for E z in terms of 
h(R g ,R) and the vertical component of the gravitational potential. 

Our next step is essentially to replace A<t> C |f by AO ac i in equa- 
tion (29) but before we do this we have to recall that the prefactor 
g(c) arose from normalising the radial probability density P(R\L Z ) 
as given by equation (6). For consistency in this formula we must 
now replace AO cff by AO a(i , with the consequence that the normal- 
ising integral is no longer analytic. Therefore our final formula for 
the distribution of at a given point in the Galaxy must be written 



n(vt\R,z) = N e - (R *- R ° }IR * 
AO ad 



2nR„K 



x exp 



tr 

f(z,R g -R), 



(36) 
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Figure 4. Demonstration of the effect on velocity distributions above the 
Sun of adding AE Z to the effective potential to allow for the tendency of 
vertical energy to shift orbits outwards. The velocity distributions are de- 
rived from fits to the models discussed in Section 4.1. Full curves show the 
full corrections with E. ("corr"), dashed lines give the same models with 
the simple adiabatic approximation ("sAA") ignoring energy conservation. 



where K(c,E,,R g ) is the numerically-determined result of normal- 
ising the revised form of P(R\L Z ) [cf eq. (12)], N normalises the 
distribution in and / is defined by equation (30). 

Fig. 4 illustrates the effect that the inclusion of AE Z in the ef- 
fective potential has on n(y$ \R, z) by showing for R = R and several 
values of z the velocity distributions predicted with (full curves) and 
without (dotted curves) AE Z . Including AE Z moves all velocity dis- 
tributions to lower v<j, particularly on the left side. The magnitude 
of the shift increases with z because at low z, n(v$\R, z) is dominated 
by orbits with small E- and therefore small AE Z . 

Far from the plane (z ~ 2 kpc) equation (36) breaks down be- 
cause the physical assumptions on which it depends fail. A problem 
that must be encountered at some height is that the vertical fre- 
quency becomes comparable to the horizontal frequency, so the as- 
sumption of adiabatic invariance of J z fails - the coupling between 
the horizontal and vertical motions becomes strong and complex. 
However, BM11 did not encounter problems with the aa below 
~ z = 2 kpc. The failure we encounter here probably arises from 
equation (47) for E z , so we below explore the effect of limiting the 
energy transfer by placing an uper limit on the value of E- to less 
than (50km s" 1 ) 2 . 



3.3 Velocity moments 

From the formulae we have in hand we can calculate a variety of 
moments. Suppose the stellar population of the disc were a super- 
position of populations that have dfs of the form 



f{E R , L., E-) oc e 



(37) 



where the dependence on E, is inspired by equation (19). Then 
since this df is a Gaussian in v R , we would have (v R ) = a 1 . To 
each value of /? s and therefore we could ascribe a fixed value 
for (v R ) = a 1 , where in general lower are connected to higher 
(v R ) = cr 2 . Hence the velocity dispersion of the entire disc could be 
obtained from the weighted average 



(v R )(R,z) = J dv 4 ,n(v^\R,z)a 2 (Rv^v c ). 
Similarly the asymmetric drift would be 
v a (R, z) = v c - J ' dvj, nfy^R, z)v # . 



(38) 



(39) 



From the work of Section 3.2 we know that E R is not strictly a con- 
stant of the motion on account of the transfer of energy between 
the radial and vertical motions, so these formulae are only approx- 
imate. We shall see that they are nevertheless useful. 

The df (37) is also a Gaussian in v ; so naively we have 
(v?) = crl for the population formed by stars of a given value of 
L-. However, when calculating (v?> for the entire disc it is essen- 
tial to bear in mind the radial variation of E z implied by adiabatic 
invariance of J z (eq. 34). In effect this variation of E z causes cr z to 
vary with radius even at fixed L z . Hence an approximate expression 
for the vertical velocity dispersion of the entire disc is 

\2/2+a 



{v 2 z )(R,z) 



I 



di> n(y f \R, z)a- z 2 (R V^Jv c ) j 



(40) 



4 APPLICATIONS 

4.1 Comparison with torus models 

To evaluate the performance of the above equations we fit the veloc- 
ity distributions of the torus model of BM1 1. The df of this model 
is 

V 



f(J r ,L z ,J z ) = U(J r ,L z )x 



2nal 



-vJ : /tri 



where 



i[- U,.l. ) 



[1 + tanh(L./L )]e 



-KJr/o-f 



(41) 



(42) 



Here Q.(E Z ) is the circular frequency for angular momentum L z , 
k(L z ) is the radial epicycle frequency and v{L z ) is its vertical coun- 
terpart and 2(2? g ) is given by equation (8). The factor 1 +tanh(L./L ) 
in equation (42) is there to effectively eliminate stars on counter- 
rotating orbits and the value of Lq is unimportant provided it is 
small compared to the angular momentum of the Sun. In equations 
(41) and (42) the functions cr-(L ; ) and cr r (L z ) control the vertical 
and radial velocity dispersions. The observed insensitivity to radius 
of the scale-heights of extragalactic discs motivates the choices 



<r r {L z ) = <r M e« {R °- R ^ 
o- z (L z ) = o- z0 ^- R ' )IR \ 



(43) 



where q = 0.45 and <t,-o and <t- are approximately equal to the 
radial and vertical velocity dispersions at the Sun. BM11 take the 
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Table 1. Parameters of the torus model's df. 

Disc Rj/kpc a>o/kms~' ov-o/kms -1 Lo/kpckms~' 

Thin 2.4 27 20 10 

Thick 2.5 48 44 10 



Table 2. Parameters of the potential employed 



Component X(Ro)/M pc 2 Rj/kpc /i/kpc i? m /kpc 



Thin 36.42 2.4 0.36 

Thick 4.05 2.4 1 

Gas 8.36 4.8 0.04 4 



Component p/M Q pc 3 q y ji ro/kpc 7", /kpc 

Bulge 0.7561 0.6 1.8 1.8 1 1.9 

Halo 1.263 0.8 -2 2.207 1.09 1000 



df of the entire disc to be the sum of a df of the form (41) for 
the thin disc, and a similar df for the thick disc, the normalisations 
being chosen so that at the Sun the surface density of thick-disc 
stars is 23 per cent of the total stellar surface density. Table 1 lists 
the parameters of each component of the df. 

BM11 took the gravitational potential to be that of Model 2 
in Dehnen & Binney (1998) modified to have thin- and thick-disc 
scale-heights of 360 pc and 1 kpc (Table 2). In this model the disc 
contributes 60 per cent of the gravitational force on the Sun, with 
dark matter contributing most of the remaining force. 

4.1.1 Distribution in v,/, 

The points in Fig. 5 show the distributions of the BM11 model 
at Rq and heights up to 2 kpc, while the curves show the fits to 
these points that we obtain from equation (36) when we approxi- 
mate AO a(i with AO cff . We take a from equation (28) and at each 
location (Ro, z) we fit the given distribution independently by ad- 
justing the parameters cr , R a and h . Table 3 gives the parameter 
values obtained from the fits and also the radial velocity dispersions 
that the fits yield through equation (38) and the true radial veloc- 
ity dispersion within the model. Notice that the parameter <x rises 



z 


CO 


Ra 


ho 


<v 2 )" 2 




x 1 





27.07 


5.30 




30.9 


33.4 


0.000095 


250 


27.56 


5.39 


229.8 


32.3 


35.2 


0.000134 


500 


29.23 


5.70 


170.4 


36.4 


40.4 


0.000265 


750 


33.30 


6.25 


190.4 


44.1 


48.4 


0.000338 


1000 


38.93 


6.86 


247.7 


53.9 


55.1 


0.000063 


1250 


40.85 


6.90 


281.3 


60.0 


58.7 


0.000079 


1500 


42.21 


7.29 


294.5 


62.4 


60.7 


0.000034 


1750 


44.45 


8.26 


284.7 


63.7 


61.8 


0.000113 


2000 


49.27 


10.90 


271.5 


64.8 


63.1 


0.000363 



Table 3. Values of the parameters for the fits shown in Fig. 5, which employ 
the simple adiabatic approximation, z denotes the distance from the plane in 
parsecs, the scale-length of radial velocity dispersion, ho the local scale- 
height. We fixed the circular speed v c = 216.25 km s _1 and disc scale-length 
i?j = 2.4 kpc to the values used by BM 1 1 . For a further comparison we give 
the rms radial velocity from eq. (38) and the value of the corresponding 
parameter o^bmii of the torus model. All fits and their x 1 values were 
derived for 60kms~' < Vd, < 260kms~'. Note that the x 1 values are not 
sensible per se, but a mere description of the relative fit quality, as there are 
no proper errors on the theoretical distributions underlying the fit. 
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Figure 5. Fitting velocity distributions of a torus model from BM1 1 at dif- 
ferent altitudes using eq. (29). Points show values from the torus model and 
the curves show our fits to these points. The scale-length of the disc was 
fixed at the BM11 value of R L \ = 2.4 kpc, the parameter a used to estimate 
the adiabatic invariant was set via equation (28). 



from 27 km s -1 at the plane to 48 km s _1 at z ~ 2 kpc and that these 
values coincide with the values of the corresponding parameters in 
the dfs of the thin and thick discs, respectively. Fig. 6 and Table 
4 show the fits obtained to the same data when AO a( j in equation 
(36) is evaluated from equation (33) with AE Z replaced by E. from 
equation (47). 

In both Figs. 5 and 6 the quality of the fits is excellent, so the 
inclusion of AE- improves the optimum fit only marginally. How- 
ever, inclusion of AE- does change the optimum value of h signif- 
icantly and in the sense of bringing it closer to the true scale-height 
of the model disc, which increases from small values very close to 
the plane, where the gas disc dominates the gravitational potential, 
through 300 pc at z ~ 300 pc, where the thin disc accounts for the 
majority of stars, to ~ 1 kpc at large heights, where the thick disc is 
dominant. Note that far from the plane the dominant population's 
scale-height h is significantly smaller than the locally measured 
scale-height of the disc because the population is dominated by 
stars withRg < R, which by equation (25) have h(R g ,R) > ho- Even 
so, when AE- is omitted, the fitted values of h are unexpectedly 
small. Including AE- increases h at all heights, while limiting the 
values of E- employed to less than (50 km s -1 ) 2 yields intermediate 
values of ho, which are not far from constant as we would wish. 
The reason adding AE- to the effective potential increases h is that 
AE Z increases the contribution to the velocity distribution at Ro of 
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ho 






27.06 


5.31 


50.5 


30.9 


33.4 


0.000097 


250 


27.46 


5.48 


347.2 


32.1 


35.2 


0.000310 


500 


29.65 


5.61 


323.4 


37.1 


40.4 


0.000267 


750 


33.99 


6.14 


326.8 


45.4 


48.4 


0.000338 


1000 


39.81 


6.79 


438.7 


55.4 


55.1 


0.000060 


1250 


41.99 


6.81 


533.2 


62.3 


58.7 


0.000085 


1500 


43.08 


7.40 


564.1 


63.6 


60.7 


0.000040 


1750 


43.83 


8.50 


496.4 


62.3 


61.8 


0.000117 


2000 


43.94 


10.70 


416.2 


58.3 


63.1 


0.000234 





27.06 


5.31 


50.5 


30.9 


33.4 


0.000097 


250 


27.46 


5.48 


347.2 


32.1 


35.2 


0.000310 


500 


29.65 


5.61 


323.4 


37.1 


40.4 


0.000269 


750 


33.99 


6.14 


326.8 


45.4 


48.4 


0.000338 


1000 


39.88 


6.74 


391.0 


55.7 


55.1 


0.000062 


1250 


41.74 


6.79 


401.8 


61.9 


58.7 


0.000081 


1500 


42.72 


7.28 


390.7 


63.4 


60.7 


0.000035 


1750 


43.69 


8.13 


345.8 


63.1 


61.8 


0.000105 


2000 


45.87 


10.18 


301.9 


61.7 


63.1 


0.000230 



Table 4. Fit parameters when the vertical energy correction AE Z is included. 
The upper half applies the unlimited correction and its fits are presented in 
Fig. 6. The lower half includes an upper limit AE. a E z = (50kms"') 2 
that produces comparable fits, but prevents a breakdown of the horizontal 
dispersion starting around z = 2kpc. Fits were taken with the same fixed 
parameters and in identical range as the fits described in Table 3, but this 
time with non-zero AE- . 



stars with small values of R g and therefore and thus reduces the 
need to suppress the contribution of the population with R s Rq, 
which dominates the peak of the distribution, relative to the stars 
that form the prominent left wing of the distribution. The other im- 
provement effected by including AE Z is to lower (v^) 1 ' 2 slightly and 
thus bring it closer to its true value at high altitudes. When there is 
no upper limit on the values of E z used in the calculation of AE Z , 
this lowering of (v 2 ,} 1 ' 2 becomes excessive above z — 2kpc because 
at such altitudes the vertical energy becomes comparable to the ra- 
dial energy. 

Irrespective of whether AE- is used, the scale-length R a ex- 
hibits a continuous rise in the fits, moving away from the value 
Rd/0.45 of the corresponding parameter of the torus model. R lT 
comes closer to R^/0.45 the lower a is chosen at higher altitudes. 



4.1.2 Distributions in Vr 

As we remarked in Section 3.3, the df (37) is such that stars of given 
L z have a Gaussian distribution in v R . Consequently, the distribution 
in Vr of all stars found at a given distance from the plane should in 
this picture be a weighted sum of Gaussian distributions with the 
weights implicit in equation (38). Fig. 7 compares this prediction 
(lines) at several altitudes z with the corresponding distributions 
from the torus models (data points). The overall agreement between 
the data points and the predictions of the formula is remarkable 
when one bears in mind that the curves have not been obtained 
by fitting to the data points. At low altitudes (red and green) the 
formula predicts a distribution that is slightly too sharply peaked 
and deficient in the wings. Around ~ 1 kpc from the plane the fit 
is near perfect. The agreement between the curve for 2 kpc and the 
data points at \v R \ < 50kms~' is much better in the lower panel 
than the upper panel, vindicating the use of the correction provided 
by AE.. 
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Figure 6. Fitting the velocity distributions of BM11 at different altitudes 
using eq. (36). The scale-length of the disc was fixed at their value of = 
2.4 kpc and the parameter a was assumed to be the function of z specified 
by eq. (28). 



4.2 Application to the Geneva-Copenhagen Survey 

We have fitted the full formula (36) at an altitude of z = 40 pc to 
the velocity distribution of the Geneva-Copenhagen Survey (GCS) 
of F and G stars (Nordstrom et al. 2004; Holmberg, Nordstrom & 
Andersen 2009). As sample we selected the full 13 520 objects that 
have measured space velocities. We adopted v c = 220kms~' and 
assumed that the Sun's velocity with respect to the Local Standard 
of Rest is 12.24km s - ' (Schonrich, Binney & Dehnen 2010), so 
232.24 km S" 1 was added to the published heliocentric veloci- 
ties. Given that the sample lies near to the mid-plane, where a = 2 
would apply, we adopted a = 1.5, and we also set the local scale- 
height to ho = 200 pc: when a and ho are allowed to vary when 
fitting to the data, they prove to be strongly correlated in the sense 
that low values of a enhance the contribution of populations with 
smaller i? g and thus favour larger values of ho for balance. How- 
ever, the differences between the residuals of the various best fits 
are not statistically significant. 

Fig. 8 shows two typical fits performed on the region 150 < 
v^/ km s~' < 250, which demonstrate how nicely and naturally the 
formula reproduces the non-Gaussianity of the azimuthal veloc- 
ity distribution. The fits are for scale-lengths R tT = 7.5 kpc (blue 
dashed line) and R a = 5 kpc (red solid line). The shorter scale- 
length provides the better fit at low v* and the worse fit to the high- 
velocity tail. The shorter length scale also yields the smaller value 
of the velocity-dispersion parameter, cr = 22.90 ± 0.45 km s~' ver- 
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Figure 7. Vr velocity distributions from the models with non-zero AE Z 
(lines) compared to the torus models (data points) at different altitudes z- 
The upper panel is for when E z is unlimited and the lower panel is for when 
T z < (50kms-') 2 . 



sus cr = 24.57 ± 0.48 km s . In the plane the core of the veloc- 
ity distribution is dominated by the youngest part of the thin disc, 
while the wings of the distribution will be dominated by the thick 
disc, and as we proceed from the core to the wings of the distri- 
bution stars of ever increasing age will grow in importance. Hence 
the true distribution in reflects the entire star-formation history 
of the Galaxy and we cannot expect to obtain a perfect fit to it by 
adjusting a single velocity-dispersion parameter, cr . Moreover, the 
statistics of the GCS catalogue to some extent reflect the complex 
selection biases involved in the catalogue's formation, and we have 
made no attempt to replicate these biases. Neglect of these com- 
plexities is presumably why our fit is less good than that obtained 
by Schonrich, Binney & Dehnen (2010). The presence of kinemati- 
cally hotter objects is confirmed by the estimated value for cr drift- 
ing to higher values when we expand the velocity interval on which 
we perform the fits. A couple of thin-disc scale-heights above the 
plane the population will be more homogeneous, being dominated 
by the thick disc, and it should be possible to obtain better fits by 
adjusting cr . 

Interestingly, on testing for a systematic shift in the formula 
recovered to 2 km s -1 the expected local standard of rest both from 
the GCS data and from the velocity distributions of the torus models 
at low altitudes. At higher altitudes the performance deteriorates 
due to the higher uncertainties. We found that for local stars AE- 
does not play a significant role, although it gives a bias of order 



Figure 8. Fitting the GCS velocity distribution with a two-parameter fit, 
using only a normalization constant and the local dispersion <tq. The lower 
panel shows the fit on a logarithmic scale to show the wings, while the 
upper panel presents the linear scale. We separately fit <tq for dispersion 
scale-length R a = 7.5 kpc (blue dashed line) and R a = 5 kpc (red solid 
line). 



1 km s -1 . In a forthcoming re-determination of the local standard of 
rest we will take this effect into account. 

Fig. 9 compares the GCS data with the v« distributions that 
follow from the fits to v<j. It is clear that the formula under-estimates 
the width of the v R distribution. In Section 4.1.2 we found that in 
the case of the torus model at small the fitted v R distribution 
was somewhat narrower than the true one even though the dis- 
tribution was very closely fitted. Correspondingly, in the case of 
the GCS data, the failure of the fit to the to adequately populate 
populate the wings of the distribution is mirrored in the fit to the v R 
distribution in Fig. 9 being least satisfactory in the wings. 



5 CONCLUSIONS 

The distribution of azimuthal velocities in the disc of a galaxy like 
ours is very skew and varies systematically with distance from the 
plane. Naturally one wants to be able to quantify such a distribution 
in an effective way. The traditional approach of fitting it with a su- 
perposition of Gaussians (e.g. Bensby et al. 2003; Ivezic et al. 2008; 
McConnachie et al. 2006) is unsatisfactory, both because there is 
no physical reasoning behind the use of a Gaussian when the distri- 
bution is not dominated by measurement error, and because when 
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Figure 9. Derived vr velocity distributions at the parameters of Fig. 8 versus 
data. 



a superposition of Gaussians is used, the parameters of the fit are 
neither unique nor physically informative. 

Our formula is based on the approximation that the vertical 
actions of stars are invariant as stars oscillate radially. We have re- 
fined this approximation by considering anew the impact that verti- 
cal motion has on the radial oscillations, which BM1 1 found to be a 
significant effect. Our treatment of this effect, being based on over- 
all energy conservation, is conceptually much sounder than that of 
BM11 and promises to play a valuable role in the interpretation 
of stellar velocities with rigorous dynamical models. However, we 
find that the power of our formula is only marginally improved by 
our more rigorous treatment of how vertical motion affects the ra- 
dial oscillations. 

Ultimately, our formula is just a fitting formula rather than 
a dynamical theory, even though we have derived it from dynam- 
ical considerations. Something the derivation highlights is how 
closely the horizontal and vertical motions of stars are intertwined, 
notwithstanding the adiabatic invariance of actions. Because both 
the vertical and horizontal random velocities of stars increase with 
age, as one moves away from the plane the mix of stars one sees 
fundamentally changes in the sense of increasing age and decreas- 
ing radius of birth. The cleanest way to model this phenomenon 
is by means of a df like those presented by Binney (2010), but a 
couple of computationally challenging steps are required to extract 
observationally testable velocity distribution such as n(v^) from a 
df: first a connection has to be established between ordinary phase- 
space coordinates and the isolating integrals upon which the df de- 



pends, and then one has to marginalise over two velocities. Evalu- 
ation of our formula is trivial by comparison. 

Our derivation makes it plain that no population of stars can 
simultaneously have scale-height and velocity dispersion that are 
both independent of radius: the sub-population formed by stars that 
have a narrow range of angular momenta must inevitably increase 
in scale-height and decrease in vertical velocity dispersion with in- 
creasing R, and if stars of larger angular momenta are added in 
to hold constant the scale-height, they will have to have an even 
smaller vertical velocity dispersion, so the vertical dispersion of 
the entire population will decline steeply outwards. Given this sit- 
uation, it is unwise to seek to define the thick disc in terms of a 
given scale-height and velocity dispersion, as some recent papers 
have done. 

We validated our fitting formula by using it to fit the distri- 
butions of components at several distances from the plane in 
a model with a well-defined df that included both thin and thick 
discs. Excellent fits were obtained. The values of the fitting param- 
eters varied slightly with the level of sophistication of the model 
employed, but were broadly in agreement with the values we would 
expect given the underlying df, especially when the most sophisti- 
cated approximations were used. This exercise implies that physi- 
cal significance can be attached to the values of parameters derived 
from fits to real data. Each fit to a distribution implies a model of 
the corresponding v R distribution. In our tests these models turned 
out to be very useful although showing a slight tendency to be too 
narrow at small |z|. 

We fitted the formula to the velocities of GCS stars and 
obtained good but not perfect fits for plausible values of the param- 
eters. The blemishes in these fits will arise from three causes: (i) 
the well known presence of pronounced clumping of stars in the 
(U, V) plane (Dehnen 1998), (ii) the need to model subtle selection 
effects in the GCS sample, and (iii) our formula is derived from an 
isothermal df and must encounter difficulty fitting data drawn from 
a system that is a superposition of systems with very disparate dy- 
namical temperatures. Hence in part the difficulties encountered in 
fitting the GCS data may reflect the importance at the extremes 
of the distribution of the thick disc and/or stellar halo. With a 
larger body of data, or data taken further from the plane, it might 
be profitable to fit the data to a sum of two or more instances of 
our formula. Our fit to the components yields a model of the v R 
components that is rather too narrow, in agreement with our work 
with the model based on a df. 

Our formula could be used to fit the line-of-sight velocity dis- 
tributions (LOSVDs) of galaxies in which individual stars are not 
resolved (e.g. Bacon et al. 2001). Since our formula has been de- 
rived on the assumption that the circular speed is independent of 
radius, it might fail to produce a satisfactory fit to the velocity distri- 
bution of stars in a galaxy with a distinctly non-flat rotation curve. 
The physical basis of our formula breaks down at distances from 
the plane of order 2 kpc, so failure to fit data for stars at higher 
altitudes might have no physical significance. 

We derived our formula by adapting to the three-dimensional 
world the planar df of Shu (1969). It would be interesting to adapt 
in a similar way the planar df of Dehnen ( 1 999), which Dehnen has 
argued is in certain respects superior to the Shu df. Unfortunately, 
the adaptation of a planar df is a non-trivial exercise on account of 
the intertwining of the radial and vertical motions mentioned above. 
Therefore in this paper we have confined ourselves to the Shu df, 
which proves to provide a very useful point of departure. 

Data for stars that lie at significant distances from the plane 
are now becoming available (Ivezic et al. 2008; Siebert et al. 201 1). 
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The measured space velocities of such stars contain significant er- 
rors arising from a combination of errors in proper-motion and dis- 
tance. It is essential to take proper account of these errors when 
inferring the true kinematics of the underlying populations. Our 
formula provides the natural way to do this: one fits the data to the 
result of folding the formula with appropriate distance and proper- 
motion errors. Schonrich, Asplund & Casagrande (2011) use this 
methodology to extract in an elegant way the information contained 
within the measured distribution of azimuthal velocities of stars 
that have quite large random velocities. We will shortly present 
similar analyses of samples that include a higher proportion of disc 
stars. 
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APPENDIX: ESTIMATING THE TYPICAL VERTICAL 
ENERGY 

In this appendix we estimate the typical value E z (R g ) of the vertical 
energy E Z (J Z ,R) of the stars that we encounter with specified an- 
gular momentum L z = Rv^, at a given location (R, z). Let cr z be the 
vertical velocity dispersion at (R, z). then 

E~ z (R g ) ^ \o- 2 z {R g ) + (X> z (z). (44) 

Regarding cr z , we have that to an excellent approximation the ver- 
tical Jeans equation reads (Binney and Tremaine 2008, eq. 4.271) 

d (p°~l) d<Jx 

— --7— -- = -p-r- (45) 
az az 

Neglecting the derivative of In a\ relative to that of In p this yields 



dz 



(46) 



where h = -(dlnp/dz) -1 is the local scale-height of the population. 
As we saw from equation (27) the local scale-height decreases to- 
wards lower guiding centre radii, if we assume all populations to 
have at their guiding centre radius a constant scale-height. Hence 
finally we adopt: 



— , dO- 

E z (R g ) = \h(R g ,R)-^- 



+ ® z (z). 



(47) 



Equation (47) involves the potential and its derivative at alti- 
tude z. We obtain these from a simple mass model with a razor thin 
gas layer that at the solar radius has surface density 12 M Q pc~ 2 and 
three exponential stellar components. In each such component the 
mass per unit surface area within distance z of the plane is 



£,(z) = 2,,o(l-e- z//, <). 



(48) 



The three components represent the thin disc, the thick disc and the 
halo with parameters 



(£1,0,22,0,23,0) = (30, 10,70)M Q pc- 
and 

(huh 2 ,h 3 ) = (300, 1000, 4000) pc. 



(49) 



(50) 



The halo contribution was chosen to match the local vertical po- 
tential from the adopted Dehnen potential, which was used for the 
torus models to which we compare our formalism - the modifica- 
tions required for a different radius or disc mass are simple. We ap- 
proximate the contribution ®,(z) to the potential from the ith com- 
ponent by assuming that the component is an infinite plane-parallel 
sheet. Then 



<b i (z) = 2nGI. w [z + h i (e- z " , <-l)]. 



(51) 



